Monte Carlo study of fermionic trions in a square lattice with harmonic confinement 
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We investigate the strong-coupling limit of a three-component Fermi mixture in an optical lattice 
with attractive interactions. In this limit bound states (trions) of the three components are formed. 
We derive an effective Hamiltonian for these composite fermions and show that it is asymptotically 
equivalent to an antiferromagnetic Ising model. By using Monte-Carlo simulations, we investigate 
the spatial arrangement of the trions and the formation of a trionic density wave (CDW), both in 
a homogeneous lattice and in the presence of an additional harmonic confinement. Depending on 
the strength of the confinement and on the temperature, we found several scenarios for the trionic 
distribution, including coexistence of disordered trions with CDW and band insulator phases. Our 
results show that, due to a proximity effect, staggered density modulations are induced in regions of 
the trap where they would not otherwise be present according to the local density approximation. 
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I. INTRODUCTION 

Ultracold atomic gases trapped in optical lattices pro- 
vide us with a new laboratory to investigate quantum 
many-body systems. They allow us not only to realize 
model Hamiltonians for electronic systems^ but also to 
investigate systems which have no direct counterpart in 
condensed matter physics. 

A shining example in this context is provided by three- 
component Fermi mixtures loaded into optical lattices, 
which have been theoretically investigated using a va- 
riety of approaches 3 ^—. The results support the overall 
idea that these mixtures give access to very different phe- 
nomena depending on the parameter regime under inves- 
tigation, ranging from Mott insulating behaviour for re- 
pulsive interaction s) 10 ' 1 1 to color-superfluidity and trionic 
phases for attractive interactions^—. Moreover, an in- 
teresting interplay between superfluidity and magnetism 
has been found to induce domain formation in globally 
balanced mixtures with SU (3) attractive interaction^, in 
striking contrast with the balanced two-component case. 
This tendency towards phase separation is quite general 
in the attractive case, being also present for asymmetric 
interactions in the strong loss regime^. In addition, it 
has been also pointed out— that multi-component Fermi 
mixtures can help to shed some light onto very complex 
phenomena connected to Quantum Chromo Dynamics. 

Testing these theoretical predictions in a laboratory is 
within the experimental capabilities of today, although 
current experiments are still performed without optical 
lattices. Indeed, mixtures of three different magnetic sub- 
levels of 6 Li [TU or 173 Yb [l4[, as well as mixtures of the 
two lowest magnetic states of 6 Li with the lowest hy- 
perfine state of 40 K [l5| already have been successfully 
trapped and cooled in current experiments. 

In this paper we focus on attractive interactions for the 
three-component mixture and more specifically on the 
strong-coupling limit, where the potential energy contri- 



bution is dominant. There is substantial evidence^^ 
that for strong-enough attraction the system undergoes 
a phase transition from a color-superfluid phase, where 
superfluid pairs coexist with unpaired fermions, to a so- 
called trionic phase, where three fermions from differ- 
ent components are bound together in new fermionic 
particles called trions. The formation of these three- 
body bound states poses new questions about their spa- 
tial arrangement and the formation of new phases in- 
volving trions as elementary objects. Several theoreti- 
cal results^^ suggest that trions tend to spontaneously 
break the translational invariance of the lattice into two 
incquivalent sublattices and give rise to staggered density 
modulations for a suitable range of parameters depend- 
ing on density, temperature, and dimensionality. Despite 
the fact that ultracold gases are charge neutral, we use 
the expression charge density wave (CDW) throughout 
the paper to identify this phase in analogy with the ter- 
minology used for electronic systems. 

The existence of a trionic CDW phase in the pres- 
ence of harmonic confinement has been investigated using 
density matrix renormalization group£ in D = 1, while 
results for the homogeneous two-dimensional case ob- 
tained within one-loop renormalization group and mean- 
field approaches^ suggest that the CDW could be the 
dominant instability at half-filling. By using dynamical 
mean- field theory on the Bcthc lattice in D = oo, it was 
shown that the superfluid phase is stable at half-filling 
against CDW in a small but finite region of coupling for 
weak attraction^. This suggests that the stability of the 
trionic CDW phase shows a marked dependence on the 
dimensionality. 

In the present work we consider a two-dimensional 
square lattice for very strong attraction both in the ho- 
mogeneus case and also with harmonic confinement. In 
order to restrict the number of parameters, we focus on a 
globally balanced system, where N a = N and N a is the 
total number of particles for the component a. We antic- 
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ipate that our results show the existence of various pos- 
sible scenarios for the spatial arrangement of trions once 
trapped and loaded into an optical lattice depending on 
the strength of the harmonic confinement and tempera- 
ture. 

The paper is organized as follows: in the next section 
we introduce a model Hamiltonian to describe a three- 
component mixture and derive an effective trionic Hamil- 
tonian for strong-coupling. Then we show that this ef- 
fective Hamiltonian can be mapped asymptotically on 
an antifcrromagnctic Ising model which we address us- 
ing Monte Carlo techniques. Details on the Monte- Carlo 
technique used are provided in Sec. IIII1 while the results 
concerning both the homogeneus case and in presence of 
the harmonic confinement are given in Sec. IIVI Finally, 
Sec. [V] concludes with a summary of the salient points 
of this paper. 



II. MODEL 

A three-component mixture of fcrmions loaded into an 
optical lattice can be suitably described by the following 
single-band Hubbard Hamiltonian: 

H = - J ^ 4,a C jM + U <r<r' n i,<r n h<7' 
{i,j),cr i,a>a' 
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Here Ci jCr (c| CT ) is the annihilation (creation) operator of 
fermions with hyperfine state a (a = 1,2,3) on the lat- 
tice site i and n^ a = ^i a c ia ^ s the fermionic number 
operator. J is the hopping amplitude between nearest 
neighbouring sites U a o' the on-site two particle in- 

teraction between fcrmions in different hyperfine states 
a and a' and fx a is the chemical potential for the species 
a. The harmonic confining potential is introduced us- 
ing the maximally packed radius r p = a^J N/ir as in Ref. 
|16| . where a is the lattice spacing, and its strength is 
parametrized by Vq. 

As outlined in the introduction, this model exhibits 
a rich variety of physical phenomena. Here, how- 
ever, we consider only the strong-coupling trionic phase 
where bound states of the three different components are 
formed^i^ and we can directly describe the system in 
terms of the composite trions. 

As we showed in Ref. ||, an effective Hamilto- 
nian for the trions can be derived by applying strong- 
coupling perturbation theory to the original Hamiltonian 
in Eq. (TTJ). By keeping only the leading and subleading 
terms, the Hamiltonian has the following form 

U e s = -Jeff / ] t]tj + V e g njnj 
(»>i> <»>j> 

El Vo 2\ T 
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where i, (t|) is the annihilation (creation) operator of a 
local trion at lattice site i and nj — t\ti is the trionic 
number operator. Since trions are color singlets, this 
Hamiltonian is analogous to a model for spinless fcrmions 
with nearest-neighbor interaction. As a consequence of 
the Pauli principle double occupancies for trions are for- 
bidden. The effective trionic hopping parameter J e fj, the 
effective interaction between trions in the nearest neigh- 
boring sites V a f[ and the chemical potential ^ e ff are given 
respectively by 
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where in the sum a, a' and a" are different from each 
other and z is the number of the nearest neighbors^. For 
the S'f7(3)-symmetric case these expressions simplif y 17 i 18 
to 
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where we assumed U < for attractive interactions. 

From Eqs igaj) and J4b}, it follows that J cS /V cS = 
J/\U\ <C 1 and therefore the effective Hamiltonian for 
trions is asymptotically dominated by the interaction 
term, while the kinetic term is subleading in the strong- 
coupling limit. In the asymptotic regime the effective 
Hamiltonian reduces to 
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where Vi = Vo(r.i/r p ) 2 . It is easy to realize that the 
Hamiltonian in Eq. ([5]) has a structure very similar to 
the one of an antifcrromagnctic Ising model in a magnetic 
field Bj 



Rising 1 ^ ^ S%Sj ^ ^ B% 



(6) 



where the parameters of the two Hamiltonians are related 
by 



1 = \v« 

Bi = - (fJ, e S - Vi) - Vtff 

s t = 2nJ - 1. 



(7a) 

(7b) 
(7c) 



3 



Correspondingly, the homogeneous system (Vq = 0) di- 
rectly maps onto an antiferromagnctic Ising model in the 
presence of a uniform magnetic field, while the effect of 
the trap is equivalent to a non-uniform magnetic field 
profile. A similar mapping on a spin model has also been 
used to investigate two-component mixtures with hop- 
ping imbalance 



.19 



III. METHOD 

In order to investigate the spatial arrangement of 
trions we perform Monte Carlo simulations of a two- 
dimensional lattice with M sites. By exploiting the map- 
ping established in the previous section the results can 
be related to the equivalent Ising spin model. 

The probability in the grand canonical ensemble for a 
specific configuration {nf} with temperature T = l/fcg/3 
is given by 
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The configuration space {nf} was sampled by using a 
Markov- chain approac h 20 ! 21 which is based on creation or 
annihilation of a single trion in a given lattice site. This 
corresponds to a single spin flip in the equivalent spin 
model. In order to have an efficient strategy for the con- 
figuration updates, we used the Metropolis algorithm^, 
where the transition probability p(a — > b) between two 
configurations is given by 



p(a — > b) = min 



1, 



p(b) 



p(a) 



(9) 



and the stationary probabilities p(a),p(b) for given con- 
figurations a and b are provided by Eq. ((5J). This allows 
the correct transition probabilities to be generated by 
computing only the energy differences between configu- 
ration a and b, connected by a single spin flip in the 
equivalent spin model. 

For the homogeneous system, which is translationally 
invariant, we have used periodic boundary conditions. 
In contrast, a system with harmonic confinement is not 
translationally invariant and periodic boundary condi- 
tions are unphysical. Therefore, in order to address tri- 
ons in a harmonic confinement, we used open boundary 
conditions, i.e. we set the occupation of the "missing" 
neighboring sites at the edges of the system equal to zero. 

Since the Hamiltonian Eq. ([5|) in the absence of har- 
monic confinement is equivalent to an antiferromagnetic 
Ising model in a uniform magnetic field, we have used 
known result o 23 i 24 i 26 to benchmark our simulations. We 
perform a careful finite-size scaling in order to con- 
firm that the thermodynamic regime has effectively been 



reached. All the results presented here are obtained from 
simulations performed on a 100 x 100 lattice. Moreover, 
in order to avoid autocorrelations in the Markov-chain 
sampling, we also perform a careful study of the equi- 
libration time t and of the number of measures we use 
to compute each observable. The excellent agreement of 
our simulation with previous results for the homogeneous 
system (see e.g. Fig. [5]) demonstrates the robustness of 
our simulations. Special care in the choice of the equili- 
bration time was required to avoid unphysical results in 
the non-homogeneus case. 

In order to quantitatively characterize the system, we 
evaluated several observables such as the local (nf ) and 
global average occupation 
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and the global CDW order parameter 
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Despite the fact that a nonzero value of C allows us to 
determine the existence of CDW order in our system, in 
the inhomogeneous case we cannot use this information 
to localize the regions where CDW order takes place. In 
this case we identify these regions by directly looking at 
the density profile (nf). Moreover, we further charac- 
terize the system by evaluating the connected density- 
density correlation function A(i,j), defined as 



A(i,i) 



((nf 



(nf))(nj-(nj))). 



(12) 



which provides us with useful information on the density 
fluctuations and the correlation length. 



IV. RESULTS 
A. Homogeneus Case 

We first consider a homogeneous lattice system with- 
out harmonic confinement. In this case, as we mentioned 
above, the Hamiltonian (Eq. ([S])) can be mapped onto an 
antiferromagnetic Ising model (Eq. (|6])) in the presence 
of a uniform magnetic field. Correspondingly, the results 
for the trionic and the Ising model can be mapped onto 
each other by the simple transformations (Eqs. (fTa |) - (jTc]) ) 
and the antiferromagnetic phase in the Ising model cor- 
responds to the CDW phase for the trionic system. 

In order to investigate the CDW order, we first study 
the CDW amplitude C as a function of temperature T 
for different values of the chemical potential, as shown in 
Fig. [TJ At T = we found the ground state of the system 
to exhibit staggered CDW order for < fj, e s < AV c s 
since C/0 (see Fig. In the homogeneous case the 
CDW phase is always characterized by a commensurate 
density (n T ) = 0.5, i.e. the density modulations in the 



4 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 
k B T/V eff 

FIG. 1: (Color online) CDW amplitude C for the homoge- 
neous trionic system as a function of the temperature T for 
different values of the chemical potential ^i e ff • 



two sublattices are always symmetric with respect to half- 
filling. Outside this range of chemical potentials, that is 
for jU e g < or /x g > 4Veff, the CDW order disappears 
and the ground state is trivially empty (n T ) = or full 
(n T ) = 1, i.e. in a band insulator phase. For /j, = 
and n = 4T4ff the CDW is degenerate with an unordered 
trionic phase, whose average density is uncommensurate 
and not fixed by the value of the chemical potential, being 
< (ht) < 0.5 for fi = and 0.5 < (n T ) < 1 for fx = 
4V c s respectively. This unordered phase has therefore 
an infinite compressibility and a very large degeneracy 
in the ground state at fixed density. Both these peculiar 
features of this phase clearly stem from having neglected 
the sublcading kinetic energy contribution, which would 
restore a finite compressibility and remove the ground 
state degeneracy leading to a metallic phase of trions^. 



With increasing temperature the CDW amplitude de- 
creases and vanishes at a critical temperature T c . The 
rounding of the transition and the non-zero magnetiza- 
tion beyond the critical temperature in Fig. [1] are due 
to finite-size effects. We summarize our results in the 
phase diagram in Fig. [5] For comparison we also show 
results for the Ising model from Ref. H^, obtained by the 
use of an extended Bethe approximation, and those from 
Ref. [U, where an analytical method based on the rela- 
tion between the singularities of the free energy and the 
zeros of the Ising pseudo-partition function on an elemen- 
tary circle were used. As one can see, we find very good 
agreement with our Monte Carlo simulations. We also 
compared our results with the Onsager solution^ in the 
absence of a magnetic field and again found very good 
agreement. 
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FIG. 2: (Color online) Monte-Carlo phase diagram for the 
two-dimensional homogeneous trionic system. The red line 
with circles marks the critical temperature T c for the tran- 
sition from CDW ordered trions to unordered trions (UT) 
within our Monte-Carlo approach. For comparison we also 
plot, after a suitable rescaling given by Eqs. (|7a,|l - (fTc| l . the 
corresponding results for the antiferromagnetic Ising model 
from Ref. [23] (blue squares), Ref. [24[ (dashed green line), 
and the Onsager solution—^ (black diamond). 



B. Trapped case 

Having studied the behavior of the homogeneus lattice 
system, we now consider the effect of superimposing a 
harmonic trapping potential. We always consider a two- 
dimensional square lattice with M = 100 x 100 sites. 

The total number of trions is fixed to N T = N = 
2500, while different values of the parameter Vo, which 
defines the strength of the harmonic confinement, and 
temperature T are investigated. 
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FIG. 3: (Color online) Global CDW amplitude C of the sys- 
tem in the presence of harmonic confinement as a function of 
temperature T for different values of trapping potential. 
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FIG. 4: (Color online) Density profiles (nf } along the cut y — 0.5 obtained within our Monte Carlo simulation (black line 
with dots) for fcsT G ft/14ff = 0.1 and different values of the harmonic confining potential Vo/Via = I (a), Vo/V e fi = 1-95 (b), 
Vo/V c s — 2.05 (c), and Vo/V e g = 2.50 (d). For comparison we also plot the LDA density profile (orange line). Below each 
panel we show a colour-coded map plot of the correlation function A(i,j — i) for i = (a;, 0.5) and j = (a;', 0.5). (see detailed 
explanation in the text). 



As a first step we calculate the global CDW order 
parameter C as a function of temperature T for sev- 
eral values of Vq (see Fig. [3j. According to our cal- 
culations, at low temperature the CDW order parameter 
C^0, and therefore there is always a finite region where 
CDW order takes place for all the setups we considered 
(Vb/Feff < 6.5). 

As already mentioned in Sec. IIII1 we then further 
characterize the regions where CDW order takes place 
by investigating the density profile (nf) and the corre- 
lation function A(i,j). In order to represent the two- 
dimensional system in our figures, we take a cross-section 
of the trap along the y = 0.5 line, which is the closest 
cut to the center of the system. 

In Fig. 0] we fix the temperature as fceT/V^ff = 0.1 



and present results for different values of the harmonic 
confining potential. For each setup we plot the density 
profile (nf ) along the cut y = 0.5 as a function of the 
x-coordinate of the lattice site i = (x,0.5). Below each 
panel we also show the corresponding correlation func- 
tion A(i,j — i) in a color-coded map as a function of the 
distance j — i = (x' — x, 0) of the two points along the 
line y = 0.5. 

For comparison with the full Monte- Carlo results for 
the inhomogeneus system, we also plot the density pro- 
files obtained within the local density approximation 
(LDA), which makes use of the Monte-Carlo results for 
the homogeneus system with an effective local chemical 
potential \ii = ^, e Q — Vi for each lattice site along the trap. 
Due to the presence of two inequivalent sublattices in the 
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FIG. 5: Typical configuration of the trapped system for 
Vo/Ves = 1.00 and kBT/V e g = 0.1, showing the presence 
of density fluctuations close to the edges of the CDW region. 



CDW phase, the LDA density profile is drawn by consid- 
ering pairs of neighbouring lattice sites corresponding to 
different sublattices in the homogeneus case. 

For weak confinement (Vo/V e g = 1 in Fig. 0^), the 
density profile clearly shows staggered density modu- 
lations in the trap center, centered around half-filling. 
These density modulations decay very fast at the bor- 
der of the CDW region, where the density decays to zero 
within few lattice sites. In this case the LDA provides 
a fairly good description of the density profile inside the 
CDW region, while deviations can be observed close to 
the edges. A snapshot of a typical configuration of the 
system for this setup is shown in Fig. [5] 

In order to better understand these deviations from 
the LDA predictions, one has to consider the correlation 
function plotted below the density profile in Fig. |H whose 
behaviour provides information about the density fluctu- 
ations and the correlation length. Indeed, for i = j, A 
provides the variance of the density distribution at the 
site i along the line y = 0.5, therefore giving account of 
the number density fluctuations in the trap. At the tem- 
perature under investigation, the thermal fluctuations in 
the bulk of the CDW regions are negligible, while sizable 
fluctuations can be observed at the borders of the CDW 
regions, as evident also in Fig. [S] 

The density fluctuations show a staggered behavior in 
correspondence to the majority and minority sublattices 
of the disappearing CDW pattern. Only close to the 
borders of the CDW domain there is a sizable spread of 
the correlation function with the distance j — i, where A 
also has a characteristic staggered behaviour that van- 
ishes within a few lattice sites. The apparent left-right 
asymmetry in the correlation plots is due to our choice 
in the lattice position with respect to the center of the 
trap, since the first lattice sites around the center are at 
x = ±0.5 and there is no central lattice site. At this 
point it is convenient to define the central average filling 




FIG. 6: Typical configuration of the trapped system for 
Vo/VeB — 2.50 and ksT/Ves = 0.1. The presence of a band 
insulating core surrounded by a CDW ring is evident, with en- 
hanced density fluctuations at the interface between the two 
regions. 



of the trap as 



T 

71 

center 



_ Si=(±0.5,±0.5)( n H 



(13) 



where n^ enter 



0.5 for the setup considered in Fig. 
With the increase of the trapping potential strength, 
shown in Fig. |4Jd, \4]p and HJi, there is a gradual suppres- 
sion of the density modulations in the trap center and the 
CDW region smoothly migrates away from the trap cen- 
ter evolving into a CDW ring for large values of the trap- 
ping potential. It is worth noting, e.g. in Fig. 0}}, that 
the local density in the central region of the trap fluc- 
tuates around an average density which is significantly 
higher than half-filling and approaches n^ nter ss 0.75 for 
the setup in Fig. [4jx For increasing trapping potential 
the density modulations eventually disappear in the cen- 
ter (see Fig. \$jp) , leading to an unordered trionic phase 
whose density increases with Vq. 

As evident in Figs.|3)D and0t, the LDA approach is un- 
able to properly describe the density profile in this regime 
of parameters, since within LDA the density modulations 
in the center disappear abruptly for much smaller val- 
ues of the confinement. The presence of wider stability 
regions for CDW ordering within the full Monte- Carlo 
calculation has to be explained through a proximity ef- 
fect: CDW order is induced in the trap center from the 
surrounding ring in the case n^ entel > 0.5, where for the 
homogeneous case this value of the density is too high 
to stabilize the CDW phase. We recall that in the ho- 
mogeneous case the CDW order is always commensurate 
in our findings and (n T ) = 0.5. This proximity effect 
in the trap can be clearly understood by looking at the 
correlation functions, which show the existence of non- 
zero correlations between different sites in the central re- 
gion of the trap; a feature that is omitted in the LDA 
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FIG. 7: (Color online) Monte-Carlo density profile (nf) for four different values of the harmonic confining potential Vo/V e g = 
1.70 (a), Vo/Ves = 1-95 (b), V /V cS = 2.05 (c), and V /V eS = 2.50 (d), and different temperatures k B T e g/V cS - 



description. A similar effect has also been observed for 
antifcrromagnetic order in a harmonic trap^S. 

Finally, for strong confinement (Vo/V e s = 2.5 in 
Fig. 0Ji), we found that the trions are in a close packed 
arrangement in the central region and therefore form a 
band insulating (BI) core in the trap center, surrounded 
by a CDW ring, as evident also in Fig. [5] The area of 
this ring decreases for increasing trapping potential, but 
as previously mentioned it never vanishes completely for 
the values of the confinement potential we considered in 
our calculations (Vo/V e g < 6.5). 

The evolution of the CDW regions at fixed values of 
Vo/V c f[ > 1.5 and for increasing temperature T is dis- 
played in Fig. [7] through the density profile. The CDW 
ordered region in the center melts at a temperature which 
slightly decreases as Vq is increased from Vq = 1.5 in 
Fig. Hi to V = 1.95 in Fig. [7p, while the CDW order in 
the surrounding ring is clearly much more robust against 
thermal fluctuations. Correspondingly the asymptotic 



value of the density in the trapping center n^ nter > 0.5 
at the melting point for the CDW core increases as a 
function of Vq. For a strongly confining trap (Vo = 2.5 
in Fig. [TJi), the CDW order is never present in the trap 
center, however, increasing the temperature causes also 
the BI core to melt. Further increase of the temperature 
has a limited effect on the remaining unordered trionic 
phase. 

Our results for the trapped inhomogeneus system are 
summarized in the phase diagram in Fig. [8] The black 
solid line in Fig. [S] always marks the critical temperature 
T c for complete disappearance of CDW order in the sys- 
tem, irrespective of the position of the CDW domain. De- 
pending on the harmonic confinement Vo/V e f{ and tem- 
perature ksT '/Veff we get four different scenarios for the 
spatial arrangement of trions: 

• A large region of CDW order in the trap center 
surrounded by unordered trions for small values of 
the confinement strength Vq and T <T C (gray area, 
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below black solid line and left from red dashed- 
dotted line in FigJBJ). Within this region we can 
further distinguish between two different cases: (i) 
For low T and small Vo, the average central filling 
"center = 0-5 (left from green dotted line) and (ii) 
for increasing T and Vb the average central filling is 
n center > 0-5 (right from green dotted line). At low 
but finite temperature the CDW core in the trap 
center disappears for Vo/V e ff = 2 and n£ QntCT = 
0.75 at the transition point. 

• For intermediate values of the confinement strength 
there is CDW order in a ring around the trap cen- 
ter, while in the central region and outside of the 
CDW ring the trions are in an unordered phase 
(green area below black solid line and between red 
dashed-dotted and blue dashed lines). 

• For large values of the confinement strength and 
low temperatures there is a BI phase of trions in 
the center, while CDW order takes place in the ring 
surrounding it (cyan area, below black solid line 
and right from blue dashed line). At the interface 
between the BI and the CDW ring as well as outside 
of the CDW ring, the trions are in an unordered 
phase. 

• The complete system is in an unordered trionic 
phase for high temperatures (white area above 
black solid line). 

We note that, as one can see in Fig. for T > 
the average central filling n^ enter . is a continuous function 
of the trapping potential and increases for increasing Vb, 
while at T = the occupation changes discontinuously 
from n^ entet = 0.5 to renter = 1- I n the latter case, the 
system stays in a CDW ground state until it becomes 
energetically favorable to move trions with large poten- 
tial energy from outside the CDW domain to the central 
region at the expense of the interaction energy. The crit- 
ical parameters for this phenomenon can be calculated 
analytically following the simple argument below. 

The energy of the outermost particle at the edge of the 
CDW phase is given by 



Eqdw = — Meff + Vo 



CDW 



— — Meff 



2Vn 



(14) 



where tcdw = V%r p is the CDW radius^, while the 
energy of this particle once moved to the trap center (for 
simplicity we assume a central site at (0, 0) is given by 



eff 



4F f 



err 



(15) 



From here it directly follows that Ecuw < Ebi for 
Vo/V e f[ < 2 and correspondingly it is energetically more 
favorable to be in the CDW phase, while for Vq/V c s > 2 
more particles will energetically prefer to move from the 
edges to the trap center and form a BI core in the center 
of the trap. 




BI in the center 
CDW in the ring 



FIG. 8: (Color online) Phase diagram of the two-dimensional 
trionic system in the presence of harmonic confinement. The 
black solid line marks the critical temperature T c for the com- 
plete disappearance of CDW order in the system. The gray 
area below this line and left of the red dashed-dotted line 
corresponds to CDW order in the trap center surrounded by 
unordered trions (UT). Within this region the dotted green 
line marks the transition between a CDW with n^ ntor = 0.5 
and n^ cntcI > 0.5. The green area below the black solid line 
and between the red dashed-dotted and blue dashed lines cor- 
responds to CDW order in a ring, while inside and outside of 
the ring the trions are in an unordered trionic phase. The un- 
ordered trionic region at the trap center evolves into a band 
insulator (BI) within the cyan area, underneath the black 
solid line and right of the blue dashed line. At the interface 
between the BI and CDW regions as well as outside of CDW 
ring the trions are in an unordered phase. Finally, above the 
black solid line the whole system is in an unordered trionic 
phase since the CDW order disappears due to thermal fluc- 
tuations. 




FIG. 9: (Color online) Average central filling n contcr as a func- 
tion of the harmonic confinement potential Vo /V e s for differ- 
ent temperatures. 
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At zero temperature and for very high trapping po- 
tential, we can expect the CDW-ring to disappear and 
the system to be in a maximally packed state with ra- 
dius r p . The energy needed to move a single particle one 
step further out from the edge of the maximally packed 
region can be easily evaluated. If this required energy is 
negative, the system will not be maximally packed in the 
ground state. The energy of a particle at the edge of the 
maximally packed system sitting at radius r p is given by 

E p = aV cff + V , (16) 

with a being the average number of occupied nearest 
neighbouring sites of the particle (a ~ 2.6 for our sys- 
tem). If this particle is moved one lattice step further 
out of the packed region, its energy is given by 

Eo = (Vo/r 2 p )(r p + a) 2 =V (l + 2a/r p + a 2 /r 2 p ) , (17) 

Therefore, the condition for the CDW-ring to completely 
disappear is given by 

AE = E p -E o <0 ^L>^_ (i 8) 

' p ' P 

At zero temperature and N = 2500 particles, r p /a = 
y / 2500/7r ?s 28.2 and the system would be maximally 
packed with no CDW-ring for a trapping potential of 

^>36, (19) 

which is about 5 times higher than the trapping poten- 
tials that we considered in our calculations. 



V. CONCLUSION 



equivalent to an antifcrromagnctic Ising model in a uni- 
form magnetic field. The results of our approach were 
found to be in very good agreement with previous re- 
sults for the Ising spin model. They show that the trions 
are arranged in a staggered density wave configuration at 
half-filling while they are unordered for an incommensu- 
rate density. This disordered phase evolves into a homo- 
geneus trionic Fermi liquid once the sublcading kinetic 
energy is reintroduced. We determined the critical tem- 
perature for the disappearance of the CDW order, which 
is analogous to the critical temperature of the equivalent 
spin model. 

In the presence of a harmonic trap we found several 
possible scenarios for the spatial arrangement of trions, 
namely coexistence of CDW domains with unordered 
trions and band insulating regions, depending on the 
strength of the confining potential and temperature. The 
CDW region moves from the center of the trap to a ring 
for increasing trapping potential. We found that the stag- 
gered density order is also induced, due to a proximity 
effect, in regions of the trap where the average density is 
not commensurate, in close analogy to the behaviour at 
the edges of the the antiferromagnetic domain in trapped 
two-component Fermi mixtures^. The inclusion of spa- 
tial correlations is crucial for a proper description of this 
feature, which is indeed missing in a LDA description of 
the system. 

An important question arising from this investigation 
is how the behaviour of the trions is modified by the 
inclusion of a finite hopping. A mapping onto an equiva- 
lent spin Hamiltonian would lead, however, to a quantum 
spin models, and therefore require an extension from a 
classical to a quantum Monte Carlo approach. For this 
reason we postpone this interesting question to a future 
work. 



In this paper we have investigated strongly attrac- 
tive three-component Fermi gases loaded into an opti- 
cal lattice and have explicitly taken the effect of a har- 
monic confinement into account. We derived an effective 
Hamiltonian for the trionic phase using a strong-coupling 
perturbation approach and showed that the trionic hop- 
ping can be safely neglected in the asymptotic regime, 
where the effective Hamiltonian describes immobile tri- 
ons with nearest-neighbor interaction. Since this model 
has a structure analogous to the antiferromagnetic Ising 
model, we performed classical Monte-Carlo simulations 
to study the spatial arrangement of the trionic particles. 

First we considered a homogeneous system without 
harmonic confinement. In this case the Hamiltonian is 
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